% Test the function makeppft filter.
%
% Yoel Shkolnisky, December 2010.

global verbose
precond=0;
precision='double';
for n=2:2:16
    A=rand(n,n);

    tic;
    [pp1 pp2]=optimizedPPFT2D(A,precision);
    B=optimizedAdjointPPFT2D(pp1,pp2,precision);
    t1 = toc;
    L=3*n-2; %n+2n-1-1
    filename=sprintf('filters%sHhatfilter%dprecond%d.dat',filesep,n,precond);
    fid=fopen(filename,'r');
    if fid>0
        N=fread(fid,1,'int');
        timestamp=fread(fid,1,'real*8');
        if verbose
            fprintf('loading 2D Hhat filter which was created on %s for n = %d and precond %d\n',datestr(timestamp),n,precond);
        end
        Hhat=fread(fid,round(L*(L+1)/2),'real*8');
        Hhat=Hhat+1i*fread(fid,round(L*(L+1)/2),'real*8');
        Hhat=reshape(Hhat,1,round(L*(L+1)/2));
        fclose(fid);
    else
        timestamp=now;
        if verbose
            fprintf('creating 2D Hhat filter on %s with n = %d and precond %d for future usage\n',datestr(timestamp),n,precond);
        end
        H=calculateH(n);
        Hhat=fftn(H,[L L]);
        fid=fopen(filename,'w');
        if fid<0
            error('failed to open file for writing, filter was not created\n');
        else
            fwrite(fid,n,'int');
            fwrite(fid,timestamp,'real*8');
            for i = 1:L
                fwrite(fid,real(Hhat(i,i:L)),'real*8');
            end
            for i = 1:L
                fwrite(fid,imag(Hhat(i,i:L)),'real*8');
            end
            fclose(fid);
            fid=fopen(filename,'r');
            if fid>0
                N=fread(fid,1,'int');
                timestamp=fread(fid,1,'real*8');
                fprintf('Filter found for n = %d, created on %s\n',N,datestr(timestamp));
                Hhat=fread(fid,round(L*(L+1)/2),'real*8');
                Hhat=Hhat+1i*fread(fid,round(L*(L+1)/2),'real*8');
                Hhat=reshape(Hhat,1,round(L*(L+1)/2));
                fclose(fid);
            end
        end
    end

    idx=[];
    jdx=1:L;
    kdx=L:-1:1;

    tic;
    Ahat=fftn(A,[L L]);

    for i = 1:L
        Htmp=[Hhat(idx) Hhat(jdx)];
        if i > 1
            idx=[i+1 idx+kdx(1:i-1)];
        else
            idx=i+1;
        end
        jdx=jdx(2:end)+L-i;
        Ahat(i,:)=Ahat(i,:).*Htmp;
    end
    Cc=ifftn(Ahat,[L L]);
    Cc=Cc(n:2*n-1,n:2*n-1);
    t2=toc;

    reportTestResult(n,(norm(B(:)-Cc(:))./norm(B(:))),t1/t2,precision);    
end

